Journal of Chemical Theory and Computation

Atomistic simulations of intrinsically disordered proteins (IDPs) are notoriously sensitive to force field inaccuracies, either regarding the protein or the water model, yielding inaccurate observables such as compactness, secondary structure propensities, or kinetics. The currently most widely used IDP force fields are Amber99sb-disp (A99disp) and Charmm36m (C36m). A99disp includes a new water model and thus optimized both, the protein and the water interactions. In contrast, C36m used the Tip3p water model and optimized only protein interactions. In many cases, C36m+Tip3p underestimates radii of gyration compared to FRET or SAXS experiments. Such overly compact structural ensembles are believed to arise from an imbalance between protein-protein, protein-water, and water-water interactions, which might be due to Tip3p inaccuracies. Here, we aim at re-balancing these interactions by combining C36m with the Optimal Point Charge (OPC) water model. Recently, this C36m+OPC combination showed improved accuracy for the disordered domain of the measles virus nucleoprotein. Here we present a systematic assessment, comparing C36m+OPC, C36m+Tip3p, C36m+Tip4p, C22*, A03ws, A99sb-ws, and A99disp for five IDPs, as well as a subset of those for five globular proteins, a set of disordered AGQ-repeat peptides, and the fast folding miniprotein CLN025. We compared extensive MD simulations (> 8.5 ms) with SAXS, NMR, circular dichroism, photo-induced electron transfer (PET), T-jump infrared spectroscopy, and X-ray crystallography measurements. We found that combining C36m with OPC improved accuracy over C36m+Tip3p for IDP ensembles without compromising its accuracy for globular proteins. While also the kinetics of the AGQ-peptides were more accurate for C36m+OPC, those of CLN025 folding were less accurate. Overall, C36m+OPC showed similar accuracy as A99sb-ws and A99disp, the latter is currently considered among the most accurate protein force fields.

Protein dynamics is fundamental to understand mechanisms. Although atomistic Molecular Dynamics (MD) remains the gold standard for predicting protein motions, its computational cost limits applications at large scales. Here, we present a coarse-grained (CG) simulation framework that combines Elastic Network Models (ENMs) with implicit-solvent Brownian Dynamics (BD) to efficiently generate protein conformational ensembles from native states. Benchmarking the method on more than 8,500 proteins against large datasets of atomistic MD trajectories and multi-state ensembles from Nucleic Magnetic Resonance (NMR), we show that short ENM-BD simulations can reproduce the residue fluctuation profiles and the dominant protein motions with remarkable accuracy. Correlations of residue fluctuations often exceeded 80%, with many proteins showing agreements above 90%, while Principal Component Analysis (PCA) revealed strong correspondences between the essential motions in the ENM-BD simulations and those observed in MD and NMR ensembles. Following CG-to-all-atom reconstruction and short energy minimization, the ENM-BD conformers also achieve high stereochemical quality, which enables their usage in downstream atomistic applications. Finally, we show that the stochastic dynamics emerging from BD trajectories closely aligns with the normal modes encoded in the underlying elastic network, highlighting that most of the intrinsic dynamics of proteins is already embedded in their native 3D topology.

Drug-target association and dissociation rates often determine in vivo efficacy more than affinity alone. The association rate constant, however, is computationally challenging to predict since the productive encounter is a rare event in a large translational and orientational search space, which lies beyond the reach of conventional atomistic simulations. We present PySTARC (Python Simulation Toolkit for Association Rate Constants), a GPU-accelerated Brownian dynamics engine for estimating bimolecular association rate constants. PySTARC is a Python reimplementation of the BrownDye engine that converges even small reaction probabilities on a single GPU, at a throughput that would otherwise require a large CPU cluster. The current framework resolves reactions between integration steps with a closed-form Brownian bridge, models the internal flexibility of the solute through a coarse-grained bead chain, distributes trajectories across multiple GPUs, checkpoints long runs, monitors convergence, and automates the entire workflow from input structures to rate estimates. PySTARC is validated against protein-ligand and protein-protein complexes spanning five orders of magnitude in the association rate constant, with most estimates within one order of magnitude of experiment.

Accurate prediction of mutation-induced changes in protein-protein affinity remains a central challenge in computational biophysics and protein engineering. Here we present a physics-based scoring method for predicting the effects of single amino acid substitutions on protein-protein and protein-peptide binding affinity. The method is based on MMPBSA energy calculations and implemented in a user-friendly, publicly available workflow (https://github.com/Hoecker-Lab/mmpbsa_scoring). To derive the scoring function, MMPBSA energy terms were fitted using a large and diverse dataset of experimentally measured binding affinity changes upon point mutation. The terms were then evaluated on independent systems. Specifically, we assessed the ability to predict mutational effects at the interface of SARS-CoV-2 spike protein with its cellular receptor, as well as computationally predicted Armadillo-repeat protein-peptide complexes. The approach performs well across diverse targets and shows promising predictive performance on both experimentally determined and computationally predicted structures. It matches or exceeds the accuracy of established physics-based protein design software and state-of-the-art deep learning-based scoring functions when applied to data outside their training domain. The method provides a computationally efficient and interpretable framework for prioritizing interface mutations in protein engineering and binding affinity optimization.

We present SOP-Multi-{pi}, an enhanced version of the two-bead-per-amino-acid resolution Self-Organized Polymer (SOP) model designed to investigate the conformational dynamics of Intrinsically Disordered Proteins (IDPs) and their phase behavior in biomolecular condensates. The original SOP-IDP and SOP-Multi models work well for chains in good solvents that predominantly adopt expanded conformations. However, they fail to reproduce the single-chain properties and phase-separation behavior of IDPs such as hnRNPA1 and its variants, which have a distinct pattern of aromatic residues. We therefore reparameterized the SOP-IDP force field to explicitly capture {pi}-{pi} and cation-{pi} interactions at the coarse-grained (CG) resolution. These interactions are weak but important for biomolecular recognition and for driving and stabilizing biomolecular condensates. Dimensionless rescaling prefactors ({lambda}) were applied exclusively to the aromatic {pi}-{pi} (Phe, Tyr, Trp) and cation-{pi} (Arg/Lys with Phe, Tyr, Trp) pairs. They were optimized iteratively against experimental radii of gyration of the hnRNPA1 low-complexity domain (A1-LCD) and its 16 aromatic and charge mutants. This restored a physically meaningful Tyr > Phe interaction hierarchy, and the Pearson correlation for Rg improved from r = 0.07 (SOP-IDP) to r = 0.83 (RMSD = 0.08 nm). The same parameters, without further tuning, gave r = 0.93 for 34 additional IDPs (24-198 residues) and reproduced Flory scaling ( {nu} ~ 0.588). We back-mapped microsecond-scale CG ensembles to all-atom representations and computed SAXS and Kratky profiles for direct comparison with experiments on IDPs and on the multidomain proteins hnRNPA1, TDP-43, G3BP1, and FUS. In direct-coexistence slab simulations, the model captures the phase behavior of A1-LCD condensates and their variants, with critical temperatures that correlate strongly with experiment (r = 0.88) but are somewhat underestimated in absolute terms. The model is implemented in LAMMPS, and a detailed user manual is provided for setting up single-chain and multi-chain condensate simulations of full-length multidomain proteins.

Simulations have emerged as a pillar in biophysics to understand behavior at the molecular scale, in particular for proteins such as intrinsically disordered proteins. These often show transient folding with long-lived states that are challenging to efficiently sample using conventional simulations. Here, we show that Hamiltonian replica-exchange at the level of hydrogen bonds can be used to systematically explore the conformational landscape of proteins. Investigating the IM30 protein, a partially alpha helical protein, under the PLUM coarse-grain force-field, we show that this approach straightforwardly allows characterization of local helix stability, which would be otherwise hard to obtain. For instance, we show that some regions of the conserved coiled-coil structure are less stable than others. By combining the different replicas with principal component analysis, we are able to dissect the conformational landscape of the IM30 protein. Taken together, this indicates that enhanced sampling remains an important part of simulations, allowing insights that cannot be easily obtained in traditional simulations.

Nucleic acid secondary structure thermodynamic models have great utility for the analysis and design of nucleic acid structures, devices, and systems, yet their complexity and piecemeal development over decades pose challenges to further innovation. Employing different functional forms for each loop type within a structure, existing nearest neighbor (NN) models necessitate numerous loop-specific dynamic programming recursions in order to calculate physical quantities of interest. Here, we derive and implement a new class of models based on transfer matrices, yielding a single functional form for the free energy of all loop types (as the trace of a matrix product) and a single dynamic programming recursion for calculating physical quantities of interest. This new class of models has additional attractive properties: 1) incorporating a subensemble of states describing proximal interactions within each loop, 2) approximating logarithmic entropic scaling for all closed loops, 3) offering the opportunity to model correlations between the subensembles of neighboring loops by generalizing transfer matrices to tensors, and 4) providing a basis for systematically adjusting the number of model parameters to enhance accuracy as new empirical data become available for parameter regression. We use literature experimental data to regress parameters for representative new transfer matrix models, observing accuracy as good as or better than that of existing NN models, yet using dramatically fewer parameters (e.g., dozens vs hundreds or thousands). For complex ensembles containing an arbitrary number of interacting nucleic acid strands totaling N nucleotides, the new single recursion leads to exact dynamic programs for calculating diverse physical quantities (including the partition function, equilibrium base-pairing probabilities, minimum free energy and suboptimal proxy structures, Boltzmann-sampled secondary structures, and gradient of the partition function) with O(N3) time complexity. The simplicity of the new partition function dynamic program enables overflow-safe partition function calculations using renormalization to achieve scalability for large complexes (e.g., 50,000 nt) as well as formulation using a global matrix multiplication primitive that facilitates GPU acceleration, yielding speedups of an order of magnitude relative to NN partition function algorithms for large complexes. The new class of transfer matrix and tensor models provides an elegant framework for future efforts to advance the accuracy, speed, and scope of computational nucleic acid analysis and design.

Ligand transport through proteins underlies many aspects of molecular function, but its detailed mechanisms remain largely inaccessible to direct observation. Molecular dynamics simulations could close this gap, but are computationally too demanding for practical purposes, including filtering stages of enzyme design pipelines. Simplified docking-based approaches might be informative, but sometimes the oversimplifications required for high throughput applications might lead to unreliable predictions. Here we present RosettaLigand PathFinder (RLPF), a fast and practical tool explicitly built to ensure conformational continuity of the ligand trajectory, protein and ligand flexibility, and robustness of every run. Built using a flexible RosettaScripts XML interface, it enables fine-grained control over all simulation parameters and adaptability for a wide range of use cases. Applied to several representative systems, RLPF reproduced established transport features and mutation effects, highlighting its ability to produce biologically relevant insights and serve as a filter in enzyme design campaigns. Wide customizability and well-defined scope of RLPF make it a useful addition to the family of tools to study the dynamics of interactions between small molecules and proteins.

A full characterisation of a molecular system requires its metastable states, their populations and their interconversion rates. Molecular Dynamics can in principle provide this characterisation but is limited by sampling and by force-field inaccuracies. Enhanced sampling addresses the first problem, but typically requires prior knowledge of reaction coordinates or metastable states and often sacrifices the kinetics, while maximum-caliber reweighting by Brotzakis, Vendruscolo and Bolhuis [Proc. Natl. Acad. Sci. U.S.A., 118(2):e2012423118, 2021] has been implemented to that matter on reweighting continuous pathways by imposing rate constants as kinetic constraints in two-state protein folding. Here, we introduce Anneal-MaxCal, which tackles both problems by imposing rate constants or time-resolved data on a Markov model learned along an annealed solute-tempering ladder spanning softer to physical interactions. MBAR assigns physical free energies to all 160 cells of a slow TICA subspace, against 79 visited at the physical coupling. The square-root approximation (SqRA), a discretised diffusion on this free-energy landscape, converts these populations into a detailed-balance rate matrix. Its prefactors, which act as a position-dependent diffusivity across each cell face, are learned from transitions at the softer couplings and extrapolated log-linearly to the physical coupling. This yields rates and a committor across the TICA space, although the physical rung contains a single barrier crossing. For alanine dipeptide, 665 ns of simulation give a relaxation time of 5430 {+/-} 996 ps against 4413 {+/-} 446 ps from 15.4 s of unbiased dynamics, whereas a constant prefactor is 2.5 times too slow. Maximum caliber then minimally perturbs this prior, with the posterior populations following from a generalised Doob transform. A synthetic time-resolved density with a tenfold slower back-transition lowers {Delta}F from 5.191 to 2.885 kBT , raises the product-basin population from 0.0055 to 0.0529 and slows relaxation 9.47-fold, leaving the reactant basin unchanged. Forward and backward rates accelerated 3 and 20-fold raise {Delta}F by 1.892 kBT and shorten relaxation from 4818 to 289 ps.

Understanding how conformational dynamics regulate protein complex assembly remains a central challenge in molecular systems biology. Here, we develop a multiscale modeling framework that integrates all-atom molecular dynamics (MD) simulations with rule-based kinetic modeling to investigate the formation of the PP1-GADD34-eIF2 complex, a critical regulator of the integrated stress response. Using umbrella sampling, we compute the potential of mean force (PMF) for key interactions, revealing strong thermodynamic driving forces for eIF2 binding ({Delta}G = -89.5 kJ/mol, KD {approx} 0.07 fM) and moderate affinity for PP1-GADD34 association ({Delta}G = -40.7 kJ/mol, KD {approx} 0.138 M). These energetics inform a PySB-based model that incorporates GADD34's folding state as a continuous variable ({varphi}), linking conformational transitions to holoenzyme assembly and activity. Simulations show that catalytic efficiency is maximized when GADD34 folding free energy is 20 kJ/mol and {varphi} = 0.9, reflecting a highly ordered state stabilized by actin and PP1 binding. Comparative analysis of energy-constrained and optimized parameter regimes reveals that folding-dependent assembly enhances both eIF2 dephosphorylation and information transmission, with channel capacities increasing from <0.002 bits to >5.5 bits. These results demonstrate how conformational regulation encodes system-level function and underscore the utility of multiscale models in bridging molecular energetics with dynamic biochemical control.

We present SMartini, an automated pipeline for deriving Martini 3 coarse-grained force-field parameters for arbitrary small molecules. Starting from a molecular structure or a SMILES string, the pipeline maps the all-atom molecule onto a coarse-grained bead representation and automatically fits both bonded and non-bonded parameters. Bonded parameters are obtained via Boltzmann inversion of atomistic molecular dynamics trajectories; successive rounds of coarse-grained simulation and distribution-matching updates then refine these estimates until the coarse-grained conformational ensemble reproduces the all-atom reference within tolerance. We validate the pipeline on a diverse set of small molecules including drug-like compounds, metabolites, and cofactors. Implemented as a modular sequence of scripts that can be orchestrated on high-performance computing clusters, SMartini makes automated parametrization of large ligand libraries practical and reduces the manual effort required for coarse-grained model development, enabling high-throughput coarse-grained simulations of protein--ligand systems within the Martini 3 ecosystem

Proteins function through dynamic conformational ensembles rather than a single static structure, yet accurately modeling these ensembles remains challenging. Here, we introduce TopoFlow, an evolutionarily conditioned flow-matching framework for protein conformational ensemble generation. TopoFlow uses evolutionary representations from multiple sequence alignments (MSAs) to encode global conformational heterogeneity. It combines this with structural variables from a variational autoencoder (VAE) to capture local flexibility. These complementary representations are then adaptively fused through a conditional modulation module to construct unified conditioning features that guide protein conformational ensemble generation. On the ATLAS molecular dynamics (MD) benchmark, TopoFlow improved the Jensen Shannon metrics for pairwise distances (JS-PwD), radius of gyration (JS-Rg), and time-lagged independent components (JS-TIC) by 10.73%, 5.78% and 7.75%, respectively, relative to BioEmu. Notably, this performance generalizes to experiment-based evaluation, where TopoFlow exhibits strong agreement with experimental observables, including chemical shifts (CS) and small-angle X-ray scattering (SAXS). Ablation analyses further indicated that the evolutionary representations and structural latent variables contributed complementary information to ensemble generation. Our results suggest that TopoFlow, through its integration of evolutionary and structural latent information, offers a promising strategy for generating structurally plausible and conformationally informative protein ensembles.

Lipid nanoparticles (LNPs) are a leading drug delivery platform for mRNA or other nucleic acid-based therapeutics. LNP platforms are typically composed of ionizable aminolipids, phospholipids, sterols, and lipids tagged with polyethylene glycol. The ionizable aminolipids are typically cationic at low pH, allowing them to encapsulate negatively charged mRNA. The pH at which the ionizable head groups change charge states, the apparent pKa, is crucial to understand the characteristics and formulation of LNP. However, this is difficult to measure as the apparent pKa differs from the pKa of the aminolipids in bulk solution due to the complex chemical and charge environment. In this study, we use molecular dynamics (MD) simulations with alchemical free energy perturbation to calculate the pKa shift as a function of lipid composition, treating LNP-like and vesicle-like environments. Our results show that the pH-dependent structure of the systems affects the trends in pKa shifts in complex ways, driven by factors including surface electrostatics and lipid internalization. Our all-atom simulation protocol allows us to efficiently estimate the apparent pKa of LNP-like models formed with novel ionizable aminolipids and estimate how it will change with composition.

Integrative structural biology seeks molecular descriptions consistent with complementary experimental measurements. A particular challenge arises for molecules whose broad conformational ensembles resist description by a few representative structures, such as natively unfolded proteins and fluid-phase biomembrane lipids. Although molecular simulations can generate such ensembles, their simultaneous agreement with a large number of experimental constraints is commonly prevented by model inaccuracies. Here we introduce autofocusing as an approach that integrates experimental information by automated finetuning of the physical model generating the ensemble. We demonstrate it for a fluid bilayer of the phospholipid POPC, for which no existing force field has simultaneously reproduced the complete set of 82 NMR C-H bond order parameters. Autofocusing brings all the 82 within experimental uncertainty while preserving agreement with independent small-angle X-ray scattering measurements and retaining or improving agreement with untargeted NMR relaxation rates and effective correlation times, thus supporting its description of membrane structure and dynamics beyond the calibration targets. The resulting simulations challenge a long-standing two-state interpretation of phosphatidylcholine headgroup structure, instead supporting a broad conformational ensemble poorly represented by crystal-derived conformers. These results establish autofocusing as a practical approach to constructing experimentally grounded conformational ensembles while retaining the computational tractability of classical molecular dynamics.

Structure-guided mutation design requires relating a proposed substitution to its molecular environment to anticipate its effects. Existing evaluations of large language models (LLMs) in protein engineering emphasize sequence-based prediction or tool-integrated workflows, leaving unclear whether LLMs themselves can infer mutation effects from structural evidence without dedicated predictors. We evaluate physical intuition-qualitative reasoning that connects structural observations with molecular interaction principles-using a benchmark that tests the same substitution in the same protein across different binding partners, with predictions based on local atomic coordinates and no dedicated mutation-effect predictors. In the tested comparisons, the best-performing LLM achieved group accuracy comparable to the evaluated energy predictors and Biomni, an LLM-driven tool workflow, while the methods differed in their prediction patterns. Most evaluated LLMs achieved higher group accuracy for larger affinity improvements. Structural interventions support selective dependence on mutation-site and partner information, while assessable structural claims largely agree with the supplied coordinates. These findings support direct LLM structural reasoning as a basis for developing rational mutation-design tools that connect molecular observations with partner-specific effects.

Interpreting single-molecule force spectroscopy (SMFS) data from optical tweezers with steered molecular dynamics (SMD) is hindered by a severe timescale mismatch: SMD trajectories span nanoseconds to microseconds, while SMFS experiments unfold over milliseconds to seconds, forcing SMD to pull six to nine orders of magnitude faster than experiment and systematically distorting rupture forces, molecular pathways, and kinetic intermediates. We introduce force-probe coupled steered molecular dynamics (fp-SMD), which explicitly propagates the force probe (the optically trapped bead) alongside the biomolecule within a single classical Hamiltonian description of the coupled molecule-linker-bead system. A Predict-Correct Trajectory Propagation (PCTP) scheme resolves the resulting computational bottleneck by decoupling the fast molecular relaxation from the slow, macroscopic bead motion, using metadynamics and transition-state theory to identify the kinetically resolved molecular intermediate reached at each macroscopic time step. On the cholesterol-{beta}-cyclodextrin ({beta}-CD) host-guest complex, fp-SMD reproduces the experimental rupture-force histogram directly at matched conditions, without extrapolation from standard SMD. fp-SMD further resolves how applied force reshapes the molecular escape pathway itself--cholesterol dissociates via a high-angle route at low force and a direct, small-angle route at high force--providing, to our knowledge, the first direct structural evidence for force-induced pathway switching within small host-guest systems believed to underlie catch-slip bond behavior. fp-SMD/PCTP thus offers a general route to simulating force-probe-coupled biomolecular mechanics at experimentally relevant timescales.

Antibiotic resistance driven by {beta}-lactamases has become a major threat to public health. Here, we propose a computational protocol for rapid estimation of {beta}-lactamase activity towards carbapenem antibiotics from sequence alone. Starting from AI-predicted acylenzyme complexes, only approximate transition state ensembles are sampled to compute catalytic electric fields, which strongly correlate with experimental activation barriers (R2 = 0.8). Compared with full quantum mechanics/molecular mechanics (QM/MM) reaction simulations, the protocol substantially reduces the amount of configurational sampling and thus computational cost. Our protocol enables efficient identification of carbapenem breakdown efficiency, further highlighting the important role of electric fields in enzyme catalysis.

Allosteric and cryptic binding sites provide unique opportunities in drug discovery for achieving selectivity, reducing toxicity and unlocking challenging therapeutical targets. Yet, cryptic pockets are particularly elusive in that they are difficult to capture with experimental methodologies, and therefore often remain undetected. Computational approaches for flexible cryptic site detection traditionally involve extensive unbiased molecular dynamics simulations, or more efficient but computationally still highly demanding biased MD alternatives that make use of enhanced sampling along low-frequency motions of the protein structure. Here we utilize a new variant of the recently introduced and orders of magnitude faster, MD-free COVIB algorithm to identify cryptic pockets efficiently in two challenging targets, the KRAS oncoprotein and the N-terminal domain of the STAT5B transcription factor. Using experimental and AI-modelled structures, respectively, we identify numerous cryptic sites, some of which have not yet been reported. The sites are prioritized and validated with our pharmacophore optimized photoaffinity-tagged fragment library via MS-based peptide mapping and NMR spectroscopy. Besides identifying new, targetable regions on the surface of these important oncotargets, we have established a protocol combining COVIB-based cryptic pocket detection and photoaffinity-based experimental peptide mapping as a robust approach for exploring unique and tractable cryptic pockets.

Cyclic peptides can engage extended protein interfaces, but sequence changes that improve developability often disrupt the selected target-bound geometry. Existing property-guided design methods lack a principled way to preserve this geometry during local multiproperty optimization. Here we present RF-PepTune, which separates structural eligibility from property prioritization. Before any property score is queried, RF-PepTune evaluates single substitutions and all position-distinct pairs induced by admitted substitutions to build a parent-specific, pair-complete admission graph. Pareto Monte Carlo tree search then optimizes predicted permeability, solubility, non-fouling, and low-haemolysis scores exclusively within this graph, and selected sequences undergo full structural re-evaluation. In a matched three-arm comparison of 334 parents across MDM2, GABARAP, and MCL1, joint property-structure success rose from 10.2% without admission to 51.5% with single-substitution admission and 77.5% with pair-complete admission; pair-complete gating improved success by 26.0 percentage points over single-substitution gating (95% CI, 20.7-31.4). Property-positive but structure-failing endpoints fell from 88.3% to 9.3%. In a separate target-balanced set, all 90 prioritized outputs improved all four model-native scores while satisfying prespecified structural criteria. This work establishes pairwise structural gating as a general strategy for constrained multi-objective molecular optimization.

Generative models of protein structures have accelerated scientific discovery. However, the leading structure-prediction models are trained primarily on crystal structures, limiting their ability to sample the conformations of intrinsically disordered proteins (IDPs). All-atom molecular dynamics (MD) simulations can probe the diversity of IDP conformational ensembles, but these simulations are slow and costly. Leveraging the speed of generative models and the accuracy of MD, we develop Zephyr, an SE(3) equivariant diffusion model trained on simulated IDP ensembles. The model generates the structure of a single protein sequence in two phases: a coarse-grained structure is first sampled and heavy atoms are subsequently denoised onto this structure. We evaluate our method on two fronts: ensemble accuracy and sampling efficiency. Across five test systems, the model reproduces key geometric and energetic metrics of the MD ensembles, while sampling at least 63 times faster.